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Abstract. Sterile neutrinos are thermalised in the early Universe via oscillations with the 
active neutrinos for certain mixing parameters. The most detailed calculation of this thermal- 
isation process involves the solution of the momentum-dependent quantum kinetic equations, 
which track the evolution of the neutrino phase space distributions. Until now the collision 
terms in the quantum kinetic equations have always been approximated using equilibrium 
distributions, but this approximation has never been checked numerically. In this work we re¬ 
visit the sterile neutrino thermalisation calculation using the full collision term, and compare 
the results with various existing approximations in the literature. We find a better agreement 
than would naively be expected, but also identify some issues with these approximations that 
have not been appreciated previously. These include an unphysical production of neutrinos 
via scattering and the importance of redistributing momentum through scattering, as well 
as details of Pauli blocking. Finally, we devise a new approximation scheme, which improves 
upon some of the shortcomings of previous schemes. 
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1 Introduction 

Light sterile neutrinos with masses around an eV have long been postulated as an explana¬ 
tion for the anomalous results seen in a number of short-baseline neutrino oscillation exper¬ 
iments (see, e.g., [1] for a review). Although intriguingly each experiment appears to point 
to a slightly different mass value, the preferred ranges of masses and active-sterile mixing 
angles are such that, within standard cosmology, the sterile neutrino state inevitably be¬ 
comes thermalised in the early Universe. The presence of an additional thermalised neutrino 
species in the eV-mass range is however in conflict with observations of the cosmic microwave 
background (CMB) anisotropies and the large-scale structure (LSS) distribution: the best 
limit from the ESA Planck mission and other astrophysical observations on the number of 
thermalised neutrino species currently stands at Aeff = 3.04±0.18 (68% C.I.) [2]. Thus, if the 
short-baseline anomalies are to be interpreted in terms of active-sterile neutrino oscillations 
despite their mutual tension, some mechanism to reconcile the sterile state with cosmology 
would be required. 

Several such mechanisms to circumvent the CMB/LSS limits on sterile neutrinos have 
been investigated in the past, some of which employ a large lepton asymmetry [3-7] or 
interactions among the sterile neutrinos [8-12] to suppress the production of sterile neutrinos 
in the early Universe. What these two mechanisms have in common is that they both delay 
the production of sterile neutrinos until at the earliest the neutrino decoupling epoch (T ~ 
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MeV), after which the active neutrino states cannot be fully repopulated even if the active- 
sterile oscillation probability should become large. This limits the total number of neutrinos 
populating the Universe, and hence avoids an unacceptably large energy density in relativistic 
particles. 

To accurately track the sterile neutrino thermalisation process—in the presence of large 
lepton asymmetry, self-interactions, or otherwise—requires that we solve the quantum kinetic 
equations (QKEs), which describe the evolution of the density matrix of the active and 
sterile neutrino phase space distributions in the presence of flavour oscillations and scattering 
(forward and non-forward). The formal expressions for all components of the QKEs up 
to order are well known [13, 14]. However, as the collision integrals describing the 
non-forward scattering are generally quite complicated and nonlinear, it is customary to 
approximate them in numerical solutions of the QKEs. Common approximations include 
neglecting Pauli blocking and feedback from the real-time neutrino phase space distributions, 
as well as ignoring the electron mass in the evaluation of the scattering matrix elements [7, 15]. 

In many ways these were reasonable approximations in their time. However, as obser¬ 
vational precision improves, it also becomes necessary to examine more closely their precise 
impact on the observables. A conservative and naive estimate of the error due to neglecting 
the Pauli blocking factors, for example, is ~ 10% [15] for active-sterile conversion before the 
neutrino decoupling epoch, which in view of Planck’s sensitivity to Aefr already appears to 
be pushing the limits of validity. As the conversion temperature approaches the decoupling 
epoch, the details of how the active neutrino states are repopulated and their oscillations 
to sterile states damped through collisions can only have an even larger impact; apart from 
thermalisation of the sterile neutrino and hence the A^eff parameter, the detailed repopu¬ 
lation mechanism affects in principle also the active neutrino energy spectra, which could 
subsequently alter the weak reaction rates during big bang nucleosynthesis (BBN). 

In this paper we include for the first time the full collision integral in the solution of 
the QKEs for active-sterile neutrino oscillations. For simplicity we restrict our attention to a 
two-neutrino model without lepton asymmetries or sterile neutrino self-interaction. We focus 
on oscillations between the electron and the sterile neutrinos described by a mass squared 
difference 6m? and vacuum mixing angle 6, although we will keep the equations general and 
bear in mind also the cases of sterile neutrino oscillations with muon and tau neutrinos. 
We assume the sterile neutrino to be heavier than the active neutrinos as this is the most 
interesting case for mass squared differences above ~ 0.1 eV^: current cosmological 
limits on the sum of the neutrino masses are in the sub-eV range [2, 16, 17], and would 
be violated if the active neutrinos were heavier than the sterile state for the values of 
interest. 

The rest of the paper is organised as follows. We begin with an introduction to the QKEs 
in section 2, where we describe the repopulation of the active neutrinos and the damping of 
flavour oscillations from collisions. We discuss various approximations found in the literature, 
and devise new, improved approximations that incorporate more of the relevant physics. In 
section 3 we first test our implementation of the full collision term for convergence, before 
proceeding to systematically compare the various approximate solutions with the full result. 
We find that although the deviations are smaller than expected for most neutrino mixing 
parameters, there are some systematic biases that can be significantly reduced using our 
new approximation scheme. We give our conclusions in section 4. Throughout the paper we 
employ a unit system in which c = h = ks = 1, and express all dimensionful quantities in 
powers of eV. 
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2 Quantum kinetic equations 


We consider oscillations between an active neutrino flavour where a = e, /r or r, and a 
sterile flavour Vg in an ensemble of neutrinos in the early universe. The density matrices p{k) 
encode the flavour content and coherence of the ensemble, and are conveniently expressed in 
terms of polarisation vectors {Po{k),P{k)), i.e., 


where fo{k) denotes a reference phase space density, which we take here to be the relativistic 
Fermi-Dirac distribution with vanishing chemical potential,^ and cr = {crx,cry,az} are the 
Pauli matrices. In this convention the active and sterile neutrino phase space distributions 
are given, respectively, by 


fuAk) = \[Po{k) + Pz{k)], 
Usik) = \[Po{k)-P,{k)], 


( 2 . 1 ) 


and the QKEs that govern their evolution can be written as [13, 15] 


P{k) = Y{k) X P{k) + - D{k)PT{k), 

fo[k) 


Po{k) 


Rqjk) 
fo{k) ■ 


( 2 . 2 ) 


Here, V(A:) = Vx{k)x + Vz{k)i contains the vacuum oscillation term as well as the matter 
potential from forward scattering, and the components are given respectively by 

Vx{k) = sin20, 

where Gj? is the Fermi constant, Mz the mass of the Z boson, rii the number density 
normalised to the equilibrium value, = 1, and = 1 + 4sec^ 0vu/ (^i^e + ^i^e) with the 
Weinberg angle 9w- Finally, Ra{k) and D{k) are respectively the repopulation and damping 
terms, which we describe in more detail in section 2.1. As we assume a vanishing lepton 
asymmetry, the same QKEs apply to both neutrinos and antineutrinos, and = np^. 


2.1 Repopulation and damping 

The repopulation and damping terms are integrals over the matrix elements for annihilation 
and elastic scattering processes. Beginning with equation (24) of [13], which includes Pauli 

^Henceforth we shall reserve /o to denote exclusively the relativistic Fermi-Dirac distribution with zero 
chemical potential, and use /eq to indicate an equilibrium distribution of a nonspecific form. 
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blocking and appears also in [18], we find these integrals to be 
Ra{k) = 2 tt j dlikidlipidYip 5E{kp\k'p') 

i 

- uAk)fpjp){i - fmm - MEp'))] 

+ Y^V^[i^a{k), j{p)Wa{k'), j{p')] [fu^{k')fj{Ep,){l - fuAk)){l - fj{Ep)) 
j 

-UAk)fj{Ep)il - UAk'mi - fj{Ep^))] , 

(2.4) 

and 

D{k) = TT J dUk'dllp'dllp 5Eikp\k'p') 

^ [fi(.Ek')fi{Epf){l - fpM) 

i 

+fpM{^-h{Ep,)){l-U{EE))\ (2.5) 

+ '^V‘^[va{k),j{p)Mk'),j{p')] [fy^{k')fj{Ep,){l - fj{Ep)) 
j 

+f,{Ep){l-f,{Ep,)){l-U^{k'))], 

where we have used the shorthand dllp = d^p/{2E)^, 6E{kp\k'p') = {Ek+Ep — E^i — Epi) is 
the ID Dirac delta function, the summation index i runs over all spectator neutrino flavours 
(i.e., where (5 ^ a) and the electron, while j runs in addition over all their antiparticles 
as well as the oscillating neutrino and antineutrino. The terms are 

VMp)MkMp'),d{k')] = {2Ef5^^\k+p,k'+p')NlNiNlNjS\MWa{p),h{k)\c{p'),d{k')), 

( 2 . 6 ) 

where iVj = y^l/2£'j,^ Ei denotes the energy of particle i, and 5|Mp(a(p), h{k)\c{p'),d{k')) is 
the squared matrix element for the forward process a{p)b{k) —>■ c{p')d{k'), summed (but not 
averaged) over initial and final spins, and symmetrised over identical particles in the initial 
and the final state. If two i/^s are present in the initial state, then SlMp must additionally be 
multiplied by 2 to account for the fact that i'a{k)va{p) —)•••• and i^a{p)^a{k) —)•... constitute 
two identical processes. 

The first part of both the repopulation and the damping integrals (2.4) and (2.5) pertains 
to annihilation processes, while the rest describes scattering processes. The repopulation 
integral (2.4) incorporates Pauli blocking in the form of additional multiplicative factors of 
the form (1 — fj) for every particle j in the final state of both the forward and reverse 
processes, and conforms with expectations. For the damping integral (2.5), however, Pauli 
blocking enters in a way that may not be entirely intuitive. Compared to the expression used 

^The prefactor Ni = s/xflEi used here differs from the definition in [13] due to the normalisation of 
the Dirac spinor. Our choice of normalisation gives the completeness relations u{pi)u{pi) = ifi + mi and 
v{pi)v{pi) = 2/i - mi. 
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by McKellar and Thomson [13], 

D{k) =7r J dUk'dHprdUp 5Eikp\k'p') 

y<Y^V'^[ua{k),Ua{p)\i{k'),i {p')] fp^ (p) + ^ [i^a {k),j {p) Wa{k'),j {p')]fj (Ep), 

* j 

(2.7) 

we find two modifications for each interaction process: one additional term (the “first term”) 
and new multiplicative factors in the second term. In terms of the evolution of the density 
matrix p, Pauli blocking enters the equation of motion as a multiplicative matrix factor of 
the form 6ij — pij (see equation (24) of [13]). The appearance of the hrst term can be traced 
to the off-diagonal part of this matrix factor, while the second term includes a factor from 
the matrix diagonal. Because the two corrections differ in sign, they cancel one another to 
some extent. 

A naive introduction of Pauli blocking into the damping integral (2.5) might lead one to 
miss the first term, which would result in an underestimation of the damping. However, as 
it turns out, the negative correction contained in the second term dominates anyhow when 
using equilibrium distributions, so that the effect of Pauli blocking is similar to what would 
naively be expected, albeit smaller. 

In appendix A we evaluate the repopulation and damping integrals (2.4) and (2.5) using 
the technique described by Hahn-Woernle, Pliimacher and Wong [19]. With this approach 
we can reduce the nine-dimensional integral to three dimensions. Of these it is possible to 
perform one integral analytically, but the remaining two must be calculated numerically. 

2.2 Approximation schemes 

The customary way to treat the collision terms is to assume most of the particles to be 
in equilibrium with the background photons. This simplifies the integrals so that solving 
the final expression is less numerically demanding. Assuming that Pauli blocking can be 
neglected and that all species follow equilibrium distributions except for the one being re¬ 
populated, the repopulation and damping terms evaluate in what we shall call the equilibrium 
approximation [4, 15, 20] to 

R,^{k) = P (^/o - ^ (Po + P.)) , (2.8) 

Eeq{k) = —P, (2-9) 

where P = CaC^pxT^, with x = fe/T, ~ 1.27, and k, 0.92 [20-22]. While this 
expression makes intuitive sense in terms of bringing the distribution towards equilibrium 
and coincides with the relaxation time approximation commonly found in the Boltzmann 
transport literature, it can also be derived from a firmer basis [13, 15]. 

An alternative approximation scheme is that introduced by Chu and Cirelli [7] (CC 
approximation), which is based on a second quantised approach [14]. Here, the combined 
collision (damping and repopulation) term is [6, 7] 

Pcou(^) =- ^— p°} — 2Gs{p — p^)Gs + {GI, p — p°} -I- 2Ga{p — p^)Ga) , 

( 2 . 10 ) 
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where Gs,a = diag ( 7 “^, O), with ( 7 !)^ = 3.06 and ( 7 ^’^)^ = 2.22 for scattering processes, and 
(7a)^ = 0.50 and ( 70 ’^)^ = 0.28 for annihilations [7], (k) ~ 3.15 T is the average momentum, 
and = /o Ih. The CC approximation was originally applied in [7] to the momentum- 
averaged quantum rate equations, and to adapt it for the momentum-dependent quantum 
kinetic equations we have had to scale equation (2.10) by a factor k/ (k) [6] to approximate 
the momentum dependence. Evaluating the matrix products, we Hnd 


Pco\\{k) 


T^— ( 4(7“)^(Paa - Pia) 

2 ^ {k) + ha?] Psa 


[( 7 ")^ + ( 7 ")^] Pas\ 
0 


which can be further recast as expressions for repopulation and damping, 
Rccik) = ^GlT^^ha? (/o - y(^’o + , 

Dccik) = [hff + haf] , 


compatible with equation (2.2). 


( 2 . 11 ) 


( 2 . 12 ) 

(2.13) 


2.2.1 Repopulation 

Observe that in the CC approximation only the annihilation processes contribute to the 
repopulation term (2.12). This approximation is reasonable for the momentum-averaged 
quantum rate equations where the expression was first used, but is not strictly correct in the 
momentum-dependent quantum kinetic equations where elastic scattering processes is also 
expected to contribute to the equilibration of individual momentum bins. The equilibrium 
approximation (2.8) on the other hand does take into account scattering processes. However, 
in doing so, it also sacrihces the principle of detailed balance in that scattering processes 
can only redistribute momentum but not repopulate a momentum bin without reference to 
the population of other bins. A better solution would be to separate the two contributions, 
but this would require a more advanced treatment than a simple relaxation towards one 
equilibrium distribution. 

Here, we develop a new repopulation approximation scheme that keeps the annihila¬ 
tion and scattering terms separate. We call this the A/S approximation. Looking hrst 
at annihilations, the full expression for annihilation into a and a is an integral over 

(/a/a ~ fuafpa)^ where we have ignored Pauli blocking factors. In the equilibrium approxi¬ 
mation, it is assumed that /j ~ /o for all particles i = a, a, Pa but Va- For a and a this is a 
good assumption. For Pa, however, we know that it overestimates if rip^ < 1. Thus, to 
compensate for the depletion of Pa, we adopt fp^ « riy^fo (remembering that np^ = riyh) in 
the annihilation part of the A/S approximation, while keeping /j = /o for i = a,a. 

For the scattering processes on the other hand, we must take care to preserve the 
neutrino number density at all times. We accomplish this in two ways, depending on whether 
the scattering process involves energy and momentum exchange between the i^a population 
and an external bath. 


1. Nonzero energy and momentum exchange. Scattering processes involving electrons, 
positrons and spectator neutrinos , all of which are assumed to be in equilibrium 
with the photons, serve to drive the iZa population to a thermal distribution with 
temperature equal to the photon temperature T. We therefore model the process using 
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a repopulation term of the relaxation form (2.8), but with the equilibrium distribution 
to which relaxes replaced by 


/. 


1 


eq, scat 


-« + !’ 


(2.14) 


where /x = is a pseudo-chemical potential.^ The ^ parameter can be determined in 
practice from the condition that the scattering contribution to must be zero, or 
equivalently, that the third (not the second) kinetic moments of /eq, scat and must 
be equal (because the collision rate is proportional to momentum). 

2. Vanishing energy and momentum exchange. Scattering amongst Va and i>a conserves 
both the energy and number densities of the active oscillating neutrinos. Energy con¬ 
servation implies that the Va population could relax to a thermal distribution with a 
temperature different from the photon temperature T, i.e.. 


/eq, Ua — 


1 


(^jT'va ~Uva IT'va -)- 1 


(2.15) 


Here, and Ty^ are determined from the combined condition that the scattering 
contributions to hy^ and Ny^ must both be zero, where Ny^ is the energy density 
normalised to the equilibrium value. The former constraint is equivalent to requiring 
equality between the third kinetic moments of /eq, y^ and fy ^, while the latter constraint 
calls for equality between the fourth moments. 

Thus, the full repopulation term in the A/S approximation can now be expressed as 


n. 


./o, 


^A/s(^) = C„,a /o - + P.)] + a 


^a,s 


~\~Ca 


jUy 


/eq, scat ~ 


-^(Po + P.) 


feq,Uc 


(2.16) 

where Ce,a = 0.180, Ce,s = 0.718, Ce,y = 0.407, = 0.102, = 0.407, and = 

0.407. Full expressions for the coefficients can be found in appendix B. As it turns out, 
the separation of the scattering processes into vanishing and non-vanishing energy exchange 
with an external bath is not crucial for the accuracy of the approximation; it has been 
included here mainly for consistency. If it were to be abandoned, one could simply truncate 
equation (2.16) at the second term, and obtain a new coefficient for the scattering term by 
adding together the numerical values of Ca,s and Ca,y given above. 

Finally, note that equation (2.16) does not accommodate a large lepton asymmetry; if 
one were present, it would be necessary to reexamine the assumptions behind the equilibrium 
approximation (2.8) and all subsequent approximations that lead to (2.16). 

2.2.2 Damping 

The damping term in equation (2.5) is affected by Pauli blocking in two ways as already 
discussed in section 2.1. As the negative correction dominates when considering equilib¬ 
rium distributions, Deq{k) in the equilibrium approximation (2.9) tends to overestimate the 


pseudo-chemical potential appears with the same sign for both particles and anti-particles, whereas the 
normal chemical potential has opposite signs for particles and anti-particles. 
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amount of damping. Conversely, Chu and Cirelli [7] included only the diagonal Pauli block¬ 
ing terms when calculating the numerical coefficients for Dccik) in equation (2.13), thereby 
underestimating the damping.^ 

Here, we propose to evaluate the damping integral (2.5) again with the approximations 
fi ~ /o for i / and ~ Due to Pauli blocking, these approximations 

lead to a damping term that is quadratic in : 

Da/S = + nu^Ca,! + Cafl) , (2-17) 

where Ce ,2 = -0.020, Ce,i = 0.569, Ce^o = 0.692, 2 = -0.020, = 0.499, and 

0 = 0.392, the negative Ca ^2 values reflecting the origin of the term in the negative 
part of the Pauli blocking factors (expressions for the coefficients are given in appendix B). 
For convenience we shall continue to call this the A/S approximation, although, unlike in the 
repopulation treatment, there is no strict separation between annihilation and scattering; the 
coefficient Ca,i is predominantly from annihilation, while comes mainly from scattering. 

Lastly, we note that it is also possible to capture some of the ^-dependence of the 
coefficients by omitting the ^-integral in the computation of the damping coefficients (see 
appendix B), and instead fit the result with a function of the form 

Cq fit = a+ he - - —, 

x + g 

where a, b, c, d, and g are constants determined by the fit. This kind of expression improves 
the agreement between the approximation and the full treatment insofar as it reproduces 
the sterile neutrino spectrum with greater success than the constant coefficients. It does 
not however lead to significant changes in AAgg and in the active spectrum which are often 
the most interesting quantities. For this reason we do not incorporate Ca,fit in the A/S 
approximation. 

2.3 Numerical implementation 

We employ a modified version of the public code lasagna^ [23] to solve the QKEs assuming 
first the full collision term (2.4) and (2.5), and then in the various approximation schemes dis¬ 
cussed above. The QKEs are solved on a fixed momentum grid, with the explicit requirement 
that P — P = 0 so as to enforce a zero lepton asymmetry. 

2.3.1 Approximation schemes 

Implementation of the approximate collision terms is straightforward in the equilibrium and 
the CC approximation. In the A/S approximation, however, additional root-finding routines 
are required to evaluate the chemical potential g of the normal scattering term (2.14), as 
well as the temperature and chemical potential /ijy^ of the r'o.t'Q.-scattering term (2.15). 
The chemical potential of the normal scattering term satisfies the equation 

gfc/T-M/T + 1^2 J fo{Po + Pz), 

■^The collision terms of the CC approximation have been presented in [7] in their integrated form, without 
details of how exactly they have been computed. However, reverse engineering suggests that they arise from the 
integral tt J dnkdUydllpidnp SE{kp\k'p')Y,iV^[i'a{k),Pa{p)\i{k'),i{p')]feq{p)feq{k){l - /eq(p'))(l “ fe<i{k')) 
for the annihilations, and tt / dUkdny dUpi dUp SEikp\k'p')J2j'^^Waik),j{p)\ua{k'),j{p')]fe<i{p)fe<i{k){l — 
/eq(A))(l “ /eq(^O) Dr the scatterings, both normalised by f dIlkfo{k). 

® Available at https://github.com/ThomasTrain/lasagna_public. 



( 2 . 18 ) 
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which can be solved numerically. The parameter space of the i^ot'o-scattering term, on the 
other hand, is two-dimensional and subject to the conditions 


/ ^ + 1 “ 2 j fo{Po + Pz), (2.19) 

/ + 1 ^2 / + -Pz). (2.20) 

In order to solve for and , we first reduce the problem to one dimension by constructing 
the ratio 

(/ dUk efoiPo + Pz))^'"" (f dHk + 1 ))"^® 

f dUk kfo(Po + Pz) ~ fdnkk/ie^/P^^-f^^^/P^^+l) ’ 


= (2-^) 


2^1/5 . 


—24Li5 ( —e 


I^Vry \ 1 4/5 


—6Li4 ( —e 


( 2 . 21 ) 


using equation (2.19) and (2.20). Here Lis(z) denotes the polylogarithm. Since the RHS 
depends on but not directly on Ty^, equation (2.21) can be solved immediately for 

l^va/Ty^ using a one-dimensional root-finding algorithm. 

Once we have established y^ua/Ty^, equation (2.20) can be evaluated explicitly for the 
temperature Ty^. We find: 


/ N 1/5 

Jdxx'^fo{Po + Pz)/2\ 

^ —24Li5 j J 


( 2 . 22 ) 


where T is the photon temperature, and again we have x = k/T. Finally, we take the ^y^ 
and Ty^ values obtained from equations (2.21) and (2.22), and further tune them by solving 
the discretised moments equations (2.19) and (2.20) using a 2D Newton’s method initialised 
with these numbers. This last step ensures that the chosen fiy^ and Ty^ values do indeed 
satisfy conservation of number and energy densities; the untuned values might not have this 
desired property because of the discretisation of momentum space. In practice, however, the 
amount of tuning is quite small. 


2.4 Pull collision term 

For the full repopulation and damping terms (2.4) and (2.5), we use the following tricks to 
simplify the calculations. Consider first the repopulation integral, which we split into three 
parts: 


Ret — Ra^e T Ra,f} T Ri 


■a,O' ? 


where the second subscript on the RHS labels the contributing processes. We use equilibrium 
distributions /eq for the electrons and positrons in the processes 

Va{k)e^{p) ^ h'a{k')e^{p'), (2.23) 

J^a{k)i>aip) ^ e~{k')e'^{p'), (2.24) 
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which should be a good assumption as these particles are tightly coupled via electromagnetic 
interactions. This assumption enables the pre-evaluation (as opposed to real-time evaluation 
while solving the QKEs) of one of the two energy integrals in each of equations (A.17), (A.18) 
and (A.20), so that the contribution of processes (2.23) and (2.24) to Ra can be expressed as 


Ra,eik) =(1 - f{k)) (y dEk>Ra,e,s,lfik') + j " f {p)) 

- f{k) Q dEk,R^^,,s,2{^ - f{k')) + J dEpf{p)R^,e,a,2^ , 


(2.25) 


with 

Ra,e,s,lik,k ) ~ y ^kd!p{Ra^s,e- T ^a,s,e+) feqiEp'){l — feq{Ep)), 
Ra,e,s,2ik, k ) ~ y ^^pi^a,s,e- T ^a,s,e+)feqiEp){l — feq{Epi)), 

Ra,e,a,2ik,P) = J dE^' Ra,a,e{'^ “ /eq(Efc/))(l - feqiEp')), 


(2.26) 


where the integration kernels Ra^xik, k',p) encode the kinematics of the interaction processes, 
and can be read off equations (A.17), (A.18) and (A.20). We note also that p' is fixed by 
energy and momentum conservation once a combination of k,k',p has been specified. 

In the limit of a massless electron, Ra,e,s,i and Ra,e,a,i of equation (2.26) scale trivially 
with temperature as oc T^. However, as the temperature drops below ~ 1 MeV, the massless 
electron approximation also becomes increasingly tenuous. Reinstating a nonzero electron 
mass significantly complicates the temperature dependence of Ra,e,s,i and Ra,e,a,i] we handle 
this by pre-evaluating equation (2.26) on a grid in {k, p, T) (or (fc, k', T)), which we interpolate 
in real time using a 3D spline when solving the QKEs. 

We use the same equilibrium approximation for the spectator neutrinos vp and antineu¬ 
trinos Pjs in the processes 


^aikivjiip) ^ Uaik'julip'), 

Va{k)h’a{p) ^ Vp{k')Vp{p'), 


where a ^ (3, andV^ are assumed to be non-oscillating. Thus Ra,p and the associated kernels 
Ra,i 3 ,s,i and Ra,0,a,i are, save for the e —>■ /3 relabelling, formally given by equations (2.25) 
and (2.26) respectively, and the corresponding expressions for Ra,xik,k',p) can be read off 
equations (A. 13), (A. 14) and (A. 19). Since the massless neutrino limit always holds in our 
considerations, the temperature dependence of Ra,0,s,i and Ra,p,a,i is trivial, and the integrals 
need only be pre-evaluated on one 2D {k,p) (or {k,E)) grid. 

Finally, for those processes involving only the active oscillating neutrinos and/or an¬ 
tineutrinos, i.e.. 


i'a{kji'i{p) ^ I'aik'ji'iip'), 

we evaluate the full double integrals (A.15) and (A.16) without approximation, since the 
momentum distributions of and Pa are a priori unknown quantities. 
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For the damping term, we see from equations (2.4) and (2.5) that it differs from the 
repopulation term only in the missing factors of f{k) and 1 — f{k). This means that the 
same simplifications of the repopulation integral discussed above and consequently all of 
the pre-evaluated integrals apply also to the damping term. For example, contribution of 
processes (2.23) and (2.24) to Da can be expressed as 


Da,e{k)=(^j dEk>Ra,e,s,lfik') + j (i£'pi2a,e,a,l(l “ /(p))^ 

- (/ dF;fc'-Ra,e,s,2(l - f{k')) + j dEpf{p)Ra,e,a,2^ , 
where Ra,e,s,i and Ra,e,a,i are exactly the pre-evaluated quantities of equation (2.26). 


(2.27) 


3 Numerical Results 

The main quantity we wish to study is the change in the energy density of neutrinos due to 
a sterile neutrino population, 

120 f 

AlVeff = - 1 = J dk{UM + USm"" ” 1- (3-1) 

This one quantity affects directly the Hubble expansion rate, making it possible to constrain 
AAgfT using various cosmological observations. 

3.1 Numerical convergence 

An important concern in the solution of the QKEs is detailed balance. If detailed balance is 
not fulfilled, at least so to a very good approximation, it will lead to unphysical excitation 
of the sterile neutrinos. As discussed in section 2.2.1, detailed balance depends on the 
implementation of and approximations applied to the repopulation term; it is violated, for 
example, by the equilibrium approximation already at the level of the equation. 

The full repopulation term, even in the presence of simplifications introduced in sec¬ 
tion 2.3, preserves detailed balance in principle. In practice, however, numerical solution of 
the QKEs using the full collision term requires that we sum over a set of discretised mo¬ 
mentum bins at every time step. This discretisation can potentially violate detailed balance, 
unless the number of momentum bins is sufficiently large. In figure 1, we solve the QKEs 
for the normalised neutrino number density n^, and effective energy density using the 
benchmark mixing parameter values 5m? = 0.1 eV^ and sin^ 29 = 0.025, employing variously 
40, 80, 100 and 150 bins; the outcomes are expressed relative to the 200 bin results. Be¬ 
yond 80 bins, convergence towards the 200 bin results to better than 0.002 across the whole 
temperature range of interest is immediately manifest. 

Comparing the different bin choices, the offset at high temperatures originates simply 
in the improved representation of the distribution function as we increase the number of bins. 
The additional discrepancy at T < 10 eV is likely to have arisen from a very small deviation 
from detailed balance, since this is temperature regime at which the thermalisation process is 
most efficient. For the remainder of the analysis we use 100 bins which gives an absolute error 
of ~ 0.001 for the benchmark mixing parameter values. This choice is a compromise between 
accuracy and speed, as the evaluation of the collision integrals scales with the number of 
momentum bins cubed: higher resolutions rapidly become prohibitively expensive in terms 
of CPU time. 
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Figure 1. Convergence test using the full collision term. We compare the difference in the neutrino 
number and energy densities between using 200 momentum bins and 40 (blue long dashed), 80 (red 
dot-dashed), 100 (green dashed) and 150 bins (cyan dotted). 


The full collision term sources from a variety of scattering and annihilation processes 
of the oscillating active neutrinos with electrons, positrons, spectator active neutrinos, and 
the oscillating active neutrinos themselves. Since we have computed their individual con¬ 
tributions explicitly, we can now also study their individual effects on the sterile neutrino 
thermalisation process. This is a sanity check, but also serves to gauge the level to which 
our implementation of repopulation conserves energy and number densities and hence fulfils 
detailed balance. To this end, we consider in figure 2 two scenarios without annihilation, and 
compute the corresponding changes in the neutrino energy and number densities relative to 
the no-oscillation case for the benchmark mixing parameter values. 

In the first scenario, we include all scattering processes (red solid and dot-dash lines) 
and hnd that the neutrino number density is conserved to an excellent degree, while the 
energy density drops below that of the no-oscillation case. This decrease can be understood 
as follows. Energy is removed from the oscillating active neutrino population through 
conversion to sterile neutrinos beginning at the low end of the momentum distribution. Some 
of these low-momentum active states are rehlled from higher-momentum states through Va 
scattering with electrons, positrons, and the spectator neutrinos. This effectively reduces 
the total energy contained in the combined active and sterile neutrino population, where the 
surplus energy has been absorbed into the background of e"*", e“, and^z/^. 

The second scenario (blue dashed and dotted lines) consists of only scattering processes 
amongst the active oscillating neutrinos themselves. Number density conservation is again 
satisfied to a good degree. Energy density is likewise approximately conserved; this is ex¬ 
pected, as the isolated Va population has no contact with the background plasma with which 
to exchange energy. Observe that the degree of non-conservation is a larger here than in the 
hrst scenario. This is because the evaluation of the Ra^a collision integrals (A. 15) and (A. 16) 
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Figure 2. Different contributions to the collision terms and their effects on the neutrino number and 
energy densities relative to the no-oscillations case. The red solid and dot-dash lines denote a scenario 
with only scattering and no annihilation. The blue dashed and dotted lines represent scattering only 
amongst the oscillating active neutrinos Va] the thick black dashed and dotted lines denote the same 
scenario, but now computed using 150 momentum bins instead of the canonical 100. 


for the VaUa processes are more sensitive to numerical inaccuracies owing to their nonlin¬ 
ear dependence on the distribution function /p^ (see section 2.3). Increasing the number of 
momentum bins from the canonical 100 to 150 bins (thick black dashed and dotted lines) 
improves the situation, although it is clear that we still do not have perfect detailed balance. 
Nonetheless, the induced error is small enough that we can conclude it is under control. 

3.2 Comparison of approximation schemes 

The different approximation schemes have different strengths and weaknesses as already 
discussed in section 2.2. On the one hand, the equilibrium approximation has the advantage 
that the assumptions involved have been analysed quite extensively [15]. It also allows the 
scattering processes to push the distributions towards the equilibrium form as they should, 
albeit at the cost of sacrihcing detailed balance. The CC approximation on the other hand 
enforces detailed balance as it prevents the same scattering processes from repopulating the 
active neutrino states. However, in doing so, the approximation also neglects possible refilling 
due to redistribution between different momentum states. With these considerations in mind, 
we expect the equilibrium approximation to overestimate the neutrino number and energy 
densities, and the CC approximation to underestimate them. 

These trends are reinforced and possibly enhanced by the behaviour of the damping term 
in the two approximation schemes. By neglecting Pauli blocking, the equilibrium approxi¬ 
mation overestimates the damping, while the CC approximation underestimates it through 
a missing positive Pauli blocking term discussed in section 2.2.2. The size of the damping 
term has direct consequences for ANefr as the production of the sterile neutrinos in our case 
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is non-resonant, where the effective production rate, 

r = ^ sin^(26'm)rcoiiision, (3.2) 

is directly proportional to the damping term FcoUision ~ D and the in-matter mixing an¬ 
gle 6m [24], Thus, the larger the damping, the higher the resulting AN^q, and vice versa. 

As shown in figure 3, the over- and underestimation of Ariiy and AAeff in the equi¬ 
librium and CC approximations, respectively, are exactly what transpires in our numerical 
solutions of the QKEs for the benchmark mixing parameter values. In contrast, the A/S 
approximation introduced in this work takes the best of both worlds, and, as is evident in 
figure 3, comes much closer to the full collision treatment than both the equilibrium and 
the CC approximations. The A/S approximation is also expected to overestimate somewhat 
the damping relative to the full treatment, since oscillations generally push the distribution 
functions below the equilibrium values. This is however a very small effect, and the fact that 
the resulting AAefr and Arii, values tend towards different sides of the full solutions suggests 
that the A/S damping term is not the cause of any major systematic bias. 

3.3 Electron mass and Pauli blocking 

Next we test the consequences of ignoring Pauli blocking, and of assuming mg = 0 in the full 
collision term. As shown in figure 3, the latter assumption makes no practical difference to 
Aril, and AN^g, since the conversion to sterile neutrinos takes place largely before electrons 
become nonrelativistic. Ignoring Pauli blocking, however, induces a ~ 0.02 absolute error as 
T drops below ~ 5 MeV, smaller than the naive expectation of ~ 10% estimated from the 
Pauli blocking terms in the relevant momentum range [15] . 

Note that there is a small subtlety when ignoring Pauli blocking: detailed balance can 
be satished for Fermi-Dirac distributions only when the Pauli blocking is taken into account. 
Otherwise, ignoring Pauli blocking generally demands that we replace Fermi-Dirac with 
Maxwell-Boltzmann statistics in order to fulfil detailed balance. We would however like to 
continue using Fermi-Dirac distributions, and therefore choose to enforce detailed balance in 
a different way. For all processes of the form a(p) -|- Va{k) ^ b{p')c{k'), we assume that 

f{Ep,)f{Ek>) - f{Ep)f{k) = /eq(Fp)(/eq(A:) - f (k)). (3.3) 

This is also the assumption behind the equilibrium approximation, and is exact if f{Epi), 
f{Ek>) and f{Ep) take on the Maxwell-Boltzmann equilibrium form. Not surprisingly we 
see in figure 3 that the equilibrium approximation solution follows the no Pauli blocking 
full solution quite well; the difference between them is due mainly to rounding errors in the 
numerical value of Ca in the equilibrium approximation immediately below equation (2.9). 

3.4 Impact on distribution functions 

It is also interesting to compare the active and sterile neutrino distribution functions in the 
different approximations. The distribution function of the electron neutrino affects directly 
the weak reaction rates in BBN [6], while the division of AN^g between and Vg have 
consequences for large-scale structure formation if 6m? is significantly larger than the solar 
and atmospheric mass splittings [11, 12]. Figure 4 shows this comparison for the benchmark 
mixing parameters. 

The first observation is that the full collision treatment gives rise to / > /o for the 
high-momentum active neutrinos at T = 10 MeV. This is not a numerical artefact, but 
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Figure 3. Comparison of the full treatment (black solid lines) with the equilibrium approximation 
(red dot-short dash), the CC approximation (blue long dashed), and the A/S approximation (purple 
dot-long dash) introduced in this work. Setting nie = 0 (cyan dotted) has almost no effect, while 
ignoring Pauli blocking in the full expression (green dashed) gives almost the same result as the 
equilibrium approximation. 


originates in both the t'^z/Q.-scattering and the annihilation processes in the presence of a 
step-like feature at low momenta, such as that produced by oscillations here. Consider first 


- 15 - 




























Relative active distribution Relative sterile distribution 



Figure 4. The active and sterile neutrino momentum distributions computed from the full treatment 
(black solid lines), the equilibrium approximation (red fit-short dash), the CC approximation (blue 
dashed), and the A/S approximation (purple dot-long dash) for the benchmark mixing parameters 
Srn? = O.leV^ and sin^(20) = 0.025 at three different temperatures. The distributions have been 
normalised to the relativistic Fermi-Dirac distribution /q. 


z^oZ/Q-scattering. Removing neutrinos from the low momentum modes causes the energy per 
neutrino to increase. This in turn triggers the number and energy conserving scatterings 
to push the distribution towards an equilibrium with a higher effective temperature, which 
automatically leads to / > /q. Annihilation processes on the other hand enhance / via a 
slightly different mechanism. For processes involving only states above the step, the equi¬ 
librium is maintained. For processes such as i^aikhigh)i^a{piow) ^ CL{k')a{p'), however, the 
reaction will be pushed to the left because there is a deficit of Va{p\ow) relative to a{k') and 
a{p'). Thus, annihilation can also overproduce z^q, at high momenta, leading to / > /q. 

The A/S approximation mimics both the annihilation and the scattering effects to an 
extent, and is the only approximation scheme that manages to reproduce / > /o at T = 
10 MeV, albeit only at very high momenta. The annihilation effect is captured by the 
riy^ factor in equation (2.16), while the separate treatment of z^qZ/q z/qZ/q accounts for the 
scattering effect. The latter constitutes the main role of the Ca,v term in the A/S repopulation 
approximation (2.16); in contrast, the Ca^s term with the pseudo-chemical potential ^ as the 
sole degree of freedom acts in the wrong direction: it becomes negative when the active sector 
is depleted, giving rise to a lower value of / for any choice of momentum. 
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Figure 5. The active momentum distributions at T = 10 MeV calculated using a collision term 
incorporating only scattering processes (top) and only annihilations {bottom). The difference between 
the full treatment (black solid) and the A/S approximation (purple long-dot dash) is much larger 
for the scatterings than for the annihilations. We have used the benchmark mixing parameters 
6m^ = O.leV^ and sin^(20) = 0.025, and the distributions have been normalised to the relativistic 
Fermi-Dirac distribution /g. 


Notwithstanding its success at capturing some degree of the / > /o phenomenon, the 
approximate treatment of scattering remains the main source of discrepancy between the 
full solution and the A/S approximation at T = 10 MeV. This is evident in figure 5 where 
we compare the distribution functions obtained assuming (i) only scattering, and (ii) only 
annihilation. The scattering-only result also demonstrates that these processes do not re¬ 
plenish the active neutrino population, yielding ///o ~ 0.85 for most momenta, while the 
annihilations are very effective at bringing ///o to 1. 

Apart from these effects, we find that the equilibrium approximation gives too large a 
final distribution for both the active and the sterile neutrinos, as expected from the oversized 
Rcq and terms. The CC approximation, on the other hand, comes surprisingly close 
to the real active neutrino distribution, but is short by ~ 5% for the sterile states. As 
repopulation does not directly affect oscillations into sterile states, we conclude that this 
offset must originate in the undersized damping term Dqc^ which as discussed in section 3.2 
affects directly the effective sterile neutrino production rate (3.2). 

3.5 Dependence on mixing parameters 

So far we have tested the various approximation schemes using the benchmark mixing pa¬ 
rameter values Sm? = 0.1 eV^ and sin^ 29 = 0.025. In reality, however, the errors caused by 
these approximations are also sensitive to the parameter values. 

Figure 6 shows AVefj at T = 0.1 MeV as a function of Sm? and sin^ 29 computed 
using the full collision term. Our results are generally quite similar to previous calculations 
using the equilibrium approximation [4, 21]. Note however that for large mixing angles and 
small mass squared differences, the contours deviate slightly from the straight lines found 
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Figure 6. ATVeg at T = 0.1 MeV using the full collision term. 


in [4, 21]. To highlight this deviation, we show in figure 7 the differences in AA^eff incurred 
respectively by the equilibrium, CC, and A/S approximations relative to the full solution in 
the (5m^,sin^ 20)-parameter space. 

As expected, the deviations always occur, independently of the approximation scheme, 
in a diagonal band corresponding to the transition region from AiVefj = 0 to AAeg = 1. Be¬ 
yond this common feature, however, the different approximations incur the largest deviations 
at different parameter values. 

The equilibrium approximation follows the result of the full calculation quite well for 
6m? values above O.OleV^, but overestimates AAefj by more than 0.1 at < 0.001 eV^. We 
can understand this deviation by looking at the conversion temperature. The temperature 
of maximal conversion is proportional to [21], so that low values of gener¬ 

ally correspond to low conversion temperatures. If the conversion temperature is sufficiently 
high (T ':$> 1 MeV), repopulation is rapid and AVeff is limited only by how fast can be 
produced through oscillations and collisions [24]. The effective production rate is given by 
equation (3.2), and production ceases as soon as the collision rate becomes too low. If most of 
the conversion occurs at low temperatures, however, collisions become too inefficient to sus¬ 
tain the population of the active sector and consequently for equation (3.2) —which assumes 
instantaneous repopulation—to hold, thereby causing the real AVefj contours to deviate from 
straight lines in relation to Sm"^ and sin^ 26 in figure 6. The equilibrium approximation errs 
in its overestimation of the repopulation rate, yielding almost straight AAefj contours even 
in the low Sm‘^, high sin^ 26 region. 

In contrast, the top right panel of figure 7 shows that the CC approximation generally 
underestimates AN^s, but works somewhat better at (5m^ < 0.01 eV^. For high (5m^ values 
the underestimation is due mainly to the undersized Dec damping term which diminishes 
the sterile neutrino production rate as already discussed in section 3.2. For low 6m? values, 
however, the agreement becomes better (a deviation of 0.02 at 6m? = 10“^, sin^ 26 = 10“^'^) 
because the deficiency of sterile neutrinos is compensated by an overproduction of active 
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Figure 7. Deviations in the ANeff{T = 0.1 MeV) values computed using various approximations 
from the full solutions. Top left: The equilibrium approximation. Top right: The CC approximation. 
Bottom: The A/S approximation. Note that the colour scale is linear, but the contour levels are not. 


neutrinos (similar to the overpopulation of the active sector discussed above in relation 
to the equilibrium approximation). This overproduction comes about as follows. Although 
technically the CC repopulation integral incorporates only annihilation processes, numerically 
the constant 2(7®)^/3.15 = 0.317 is rather larger than its annihilation counterpart in the A/S 
approximation, Ce,a = 0.177. Thus, the CC approximation does inadvertently contain some 
degree of repopulation due to scattering, which drives up the active neutrino repopulation 
rate relative to the true rate. 

Finally, the bottom panel of hgure 7 shows the A/S approximation. Here, the agreement 
is generally much better than either the equilibrium or the CC approximation, although 
there remains a small discrepancy of ~ 0.01 in the low dm?, high sin^ 20 corner. This is 
not surprising, as the region is characterised by large spectral distortions from conversion 
at low temperatures, and is thus most sensitive to how exactly we handle repopulation in 
the active sector. It is nonetheless remarkable that for most of the parameter region, the 
A/S approximation is able to reproduce the full results at a precision of ~ 0.001 through 
a fairly simplistic description of the collision term. This deviation is in fact comparable to 
the numerical error expected of the full treatment due to our choice of momentum resolution 
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(see section 3.1). Therefore, figure 7 not only validates the A/S approximation, but through 
comparison with a physically intuitive modelling of collisions, also confirms that the full 
collision term has been implemented correctly 

4 Conclusions 

We have calculated in this work the full collision term for active-sterile neutrino oscillations 
in the early universe, and for the specific case of (r'e, i^s)-oscillations implemented it in the 
computation of AA"efr from sterile neutrino thermalisation. In particular we have included for 
the first time a nonzero electron mass and Pauli blocking in the collision integrals. The former 
turns out to have a negligible impact on the final results, while ignoring the latter gives rise 
to noticeable discrepancies. Nonetheless, our full treatment confirms previous analyses based 
on approximations of the collision terms that a 1 eV sterile neutrino coupled with a mixing 
angle sin^ 29 ~ 0.1 produces AAefr ~ 1; the discrepancies arising from approximations are at 
the level that affects only precision calculations. 

We then proceeded to perform a systematic comparison of the full collision treatment 
with two approximation schemes found in the literature. We find that the commonly used 
“equilibrium approximation” [15] reproduces AAgff at better than 0.04 for dm? > 0.01 eV^, 
but incurs large errors (> 0.1) for very small mass squared differences due to the unphysical 
repopulation of the active sector by elastic scattering. The “CC approximation” of Chu 
and Cirelli [7] is discrepant up to 0.04 for 5m? > 0.0001 eV^, but improves for low 5m? 
values because of a fortuitous cancellation between an underestimated damping term and an 
overestimation of repopulation from annihilation processes. 

Recognising that the equilibrium and CC approximations neglect different physical ef¬ 
fects, we have devised a new approximation scheme, the A/S approximation, in which scatter¬ 
ing and annihilation contributions to repopulation are treated separately, and the damping 
term includes Pauli blocking. As the scheme is better able to capture the physics of re¬ 
population and damping, it also pushes the error in AAgfr down to 0.001 for most of the 
parameter region, although minor deviations (~ 0.01) remain in the 5m? < 0.001 eV^ region, 
where the sterile neutrino conversion temperature is lowest and hence spectral distortions 
from oscillations are expected to have the largest effect. 

The connection between low 5m? deviations and low temperature spectral distortions 
also has implications for an inverted mass spectrum, i.e., where the active state is heavier 
than the sterile state, as well as for the various mechanisms designed to reconcile eV-mass 
sterile neutrinos with cosmology. In the case of an inverted mass spectrum, sterile neu¬ 
trino thermalisation proceeds via a resonance, which, depending on the adiabaticity of the 
resonance and hence the mixing angle, can cause more disturbance to the active neutrino 
spectrum than in the non-resonant case. Likewise, mechanisms that block the production of 
eV-mass sterile neutrinos typically work by delaying the thermalisation to low temperatures, 
which by construction also makes them more sensitive to the approximations employed for 
the collision terms. 

At the current level of observational precision—cT(AAeg) ~ 0.2 from Planck [2]— our 
analysis shows that the CC and the A/S approximations can be reliably applied to most 
active-sterile oscillation scenarios, whereas the equilibrium approximation appears to be 
approaching its boundary of validity if the sterile neutrino conversion temperature is too low. 
In the future, large-volume galaxy surveys are expected to improve the sensitivity to AAgff to 
~ 0.03 [25]. Thus, should hints for a 1 eV sterile neutrino persist in the laboratory, a collision 
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treatment more precise than either the equilibrium or the CC approximation can offer would 
become necessary. Short of evaluating the full collision terms, which is numerically costly or 
possibly even infeasible in some cases, the A/S approximation developed in this work appears 
to be a most convenient alternative. 
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A Derivation of the full collision terms 

A.l Neutrino—electron scattering in the s-channel 


Consider first the case of t'o + e ^ z/q, + e . The repopulation term due to this process 
takes on the form 

j d’k'd’p'd’p «(“)((: + p-k'- ?>') 

-fM)fe{Ep){l - - fe{Ep,))] , 

(A.l) 

where the matrix element is [26, 27] 


|A/|2 = 32G| {Al{k-p){k'-p') + Bl{k-p'){k'-p) - mlAo,Ba{k ■ k')) , (A.2) 

with Ag = 2sm^6w + Ij = 2sin^0vu — B^ = 2s\v?9w-, and B^^r = 2sin^0iv Note 
the change of notation here in the appendix: p and p now denote respectively the 4- and 
3-momentum, while in the main text we use p to indicate the magnitude of the 3-momentum. 

The matrix element has three different dependences on the momentum, and therefore re¬ 
quires three different parameterisations to simplify the integral. We label the first dependence 
{k ■ p){k' ■ p') the s-channel, the second {k ■ p'){k' ■ p) the n-channel, and the last dependence 
{k ■ k') the t-channel, where for each channel it is convenient to define a 3-momentum variable, 
given respectively by 


q = p -I- k = p' -I- k', 

V = p - k' = p' - k, 

w = k — k^ = p^ — p. 


which replaces one of the integration variables of equation (A.l). Then, to evaulate the inte¬ 
gral for each channel we simply follow the technique described by Hahn-Woernle, Pliimacher 
and Wong in [19]. 

Taking the s-channel as a worked example, we use the 3-momentum variable q as a 
reference and define around it a coordinate system 


q = |q|(o,o, 1), 

k = Efc (0, sin rj, cos rj ), 
k^ = Efc/(cos())sin0,sin())sin0,cos0). 
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With this choice, we find the quantities 


s — (p + — {p' + k')'^ — {Ep + Ek)'^ — |qp, 


q-kf 

|q-k|2 



qp + |k'p — 2q-k' = |qp + E^, — 2|q|£'fc/ cos0, 
qp + |kp - 2q - k = |q|^ + EI - 2|q|Efc cosry, 


and consequently 

|M,|2 = 8 GlAl{{Ep + Ekf - |qp - mlf (A.3) 

for the first term of the matrix element (A.2). 

Since the matrix element (A.3) now depends only on energies and the magnitude of q, 
we can use the Dirac delta functions in the integral (A.l) to integrate out the directional 
dependencies. Evaluating first the p'-integral as [19] 


-p') = I d^p'dEp, ^^^^' ^ .^^9{Ep, - rue) 

^Epi J 2y'|p'|^+m| 

X 6{Ek + Ep- Ey - Ep,)5\\^ + p - k' - pO 


5{Ek +Ep- Ek, - V|q - kf + 


m; 


2y^|q-k'|2 + 


mi 


'd{Ek + Ep — Ej^i 


we use the property 




E 

f{xi)=0 


djx - Xj) 


me), 

(A.4) 


and hence 


- |piP - mj) = 


5{Ei - -y/|pi|2 + mj) 5{Ei + y/|pi|2 + m|) 


‘^\\Vi\ +mf 


2y |Pip + m‘f 

to further simplify equation (A.4) to [19] 

/ ^ P ~ ~ 

= 5{{Ek + Ep- Ek')"^ - |q - k'l^ - ml) 6{Ek + Ep - Ek’ - me) 

= 6 ((Efc + Ep — Ek’)"^ — [q]^ — E^i + 2|q|Efc/ cos0 — uig) 0[Ek + Ep — Ek’ — 
Ef.! — {Ek + Ep — Ek’)"^ + |qp + ml 


me 


6 cos 9 — 


2\q\Ek’ 

In a similar way, we rewrite 
d^p 


2|q|Efc/ 


6{Ek + Ep- Ek' - me). 

(A.5) 


2Er, 


= J d^p dEp S{Ep — |q — kp — ml) 9{Ep — me) 

El-E^^ + H-^ + ml' 


/ 


= / d qdEp 


1 


(A.6) 


2|q|Efc 


6 cos p — 


2|q|Efc 


9{Ep - me), 
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where in the second line we have also changed the integration variable from p to q. Then, 
applying equation (A.5) and (A.6) to the integral (A.l), we find, after performing the trivial 
angular integrations and averaging over the direction of the incoming neutrinos (f dcos?]/2), 
the s-channel contribution 


Ra,s,e,s ~ c£i 3^2 J d COS T] d COS 6 d\q\ dEk> dEp 6{cos9 — .. .)5{cosr] — ...) 


(A.7) 


X \Ms\‘^F9{Ek + Ep- Ek> - me)9{Ep - me), 


where E denotes all the distribution functions, and the arguments of the two Dirac delta 
functions can be read off equations (A.5) and (A.6) respectively. 

The integrals over cos 9 and cos r] involve only Dirac delta functions. Taking the cos rj- 
integral as an example and noting that integration limits can be equivalently expressed as 
step functions, we find 


J d{cos r])6 ( cos rj 


= 9 1- 


El-E^p + \ci\^ + ml ' 
2|q|£^fc 

El-E^p + \ci\^+ml'\ 


El 


2|q|.Efc 


El + |qP + ml ^ 
2|q|T’fc 


The two step functions can be reinterpreted as limits on |q|, and together they confine |q| to 
the interval 

||k| - IpII = \Ek - < |q| < Ek + El- ml = |k| + |p|. 

Integrating over d cos 9 gives the same formal result save for the replacements k — )• k' and 
p —7- p'. Thus the combined integration limits on |q| can be written equivalently as 

max (||k| - IpII, ||k'| - |p'||) < |q| < min (|k| + |p|, |k'| + |p'|) , (A.8) 

and we note that |p'| is determined from |p|, |k|, |k'| by imposing energy conservation. Ap¬ 
plying the limits (A.8) to the integral (A.7) then yields 


Re 


2 roo rc 


dEr, 


d\q.\m^E 

max(me,fc'—Ai+me) (A.9) 

q| - max(||k| - |p||, ||k'| - |p'||))6»(min(|k| |p|, |k'| |p'|) - |q|) 


as our reduced repopulation integral from the s-channel. The u- and t-channel integral 
reduction proceeds in a similar manner, using v and w respectively as an integration variable. 


A.2 The massive case 

The reduced integral (A.9) is but one of three contributions to the repopulation of t'o arising 
from neutrino scattering with electrons. The full collision term, including scattering and 
annihilation processes with has in total 14 such terms to be evaluated (see table 1). 

Fortunately, however, these 14 terms can all be recast into one of the standard s-, t-, and 
tt-forms, and thus can be handled in ways similar to that discussed above. 

As the reduction procedure concerns only kinematics and does not involve the actual 
matrix element besides the initial classihcation of the momentum dependence into s-, u- or 
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Reaction (a ^ /3) S\M\'^ 


l'a{k)l'p{p) 

-)> 

I^a{k')v 0 {p') 

i 2 GUk-p){k'-p') 


Voi{k)h'p{p) 


Va{k')vp{p') 

32GUk-p')ik'-p) 


Va{k)Va 

'.(P) 


Va{k')Va 

ip') 

32Gj,2{k-p){E-p') 


Va{k)Da 

'{P) 

-)■ 

^a{k')Da 

(P') 

32Gj,A{k-p'){k'-p) 


Va{k)e~ 

'{P) 

-)> 

Va{k')e~ 

■(/) 

32G|((2xw±1)^(A' 

■p){k' -p') + 4:x'^{k ■ p'){k' ■ p) 

— { 2 xw ± l) 2 xw'm‘l{k ■ k')) 

Va{k)e'^ 

'{P) 


Va{k')e'^ 

■{p') 

32G|,((2a:w±l)^(fc' 

■ P')ik' ■ p) + {k-p){k' -p') 

— { 2 xw ± l)2xwmg(fc • fc')) 

Va{k)Va 

dj>) 


V 0 {k')vp{p') 

32Gl{k-p'){k'-p) 

Va{k)Da 

{P) 

->• 

e {k')e'^ 

-ip') 

32Gl{{2xw±l)^{k- 

■k'){p-p') + 4,x'^{k ■ p'){k' ■ p) 

— { 2 xw ± V) 2 xw'm 1 {k ■ p)) 


Table 1. Matrix elements for all relevant reactions involving 1 /^, with xw = sin^ 6 w = 0.23864 [28], 
S' is a symmetrisation factor of 1/2 for each pair of indistinguishable particles in the initial and the 
final state, and |Mp has been summed but not averaged over initial and final spins. For the process 
with two t'aS in the initial state, we have further multiplied the matrix element by 2 to account for 
the fact that i'a{k)va{p) —>■... and Va{p)i^a{k) —?►... constitute two identical processes. Where there 
is a choice of ±, the plus signs are for a = e, and the minus signs for a = p,T. The corresponding 
matrix elements for Pq, can be obtained by the exchange {k ■p){k' -p') O {k ■p'){k' -p) for the elastic 
scattering processes. 


t-forms, we shall keep the calculation as general as possible and allow for the possibility that 
all initial and final states are massive. Then, the number of independent reductions to be 
performed is only three, which yield: 




1 


Q, s—channel — 


277r3£;.|k| 


COO coo 

/ dEk' / 

Jmu J m, 


dE^ / d\q\ S\Ms\‘^F 


t— 


f my ./max(mp,A:'—fc+rrip/) J 

X e{\q\ - max(||k| - |p||, ||k'| - |p'||))6l(min(|k| + |p|, jk'j + jp'j) - jqj), 

(A.IO) 
Ak 


Q.i—channel — 


2'^TT^Ek\k\ 


COO CC 

/ dEk' / dEp 

Jmy v/max(mp,/c'—fc+m /) 


J djwj S\Mi 




X 0(|w| - max(||p| - jp'll, ||k| - |k'||))6'(min(|p| + jp'j, jkj + jk'j) - jwj), 

(A.ll) 


^a,n-channel -27^3^^|k| 


COO cc 

/ dEk' / 

Jmy J m-. 


dE, 


J d\v\ S\Mu 


' max(mp —k-\-m^i ) 

X 0(|v| - max(||k| - |p'||, ||k'| - |p||))0(min(|k| + jp'j, jk'j + jpj) - jvj). 

(A.12) 

After inserting the matrix element S\Mx\^ and momentum distributions E, these reduced 
integrals are valid for any 2 —)• 2 process. 


A.3 The full collision terms 

The matrix elements for all elastic and inelastic processes involving Va at temperatures T < 
rrip are summarised in table 1. These have been computed at various times by several different 
groups [26, 27, 29], but can be easily obtained from first principles in the four-fermion limit. 
Using these matrix elements and the expressions (A.IO), (A.ll) and (A.12), we can now 
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determine the contribution of each process to the repopulation integral. The results are as 
follows. 


1. Va{k)vi3{p) Va{k')yp{p')-. 

poo pc 

/ dEy / 

^0 J mi 


G‘^ 

^ 2(27r)3E2 


dE, 


■L 




msociQ.Ey—Ek) JmQyi{\E]z—Ep\,\2Ey-Ep—E]z\) 


(i|q| 


X 


(Ep + ~ |q|' 


f{Ek>)f{Ep,)il - f{Ek))il - f{Ep)) - f{Ek)f{Ep){l - f{Ek^)){l - f{Ep,)) . 

(A.13) 


2. l ' a { k ) vp { p ) l ^ aik ') l '^ ip '). 
Gl 




2(27r)3i?2 


poo p 

/ dEk' 

Jo Jrc 


dEr. 


inciiu.(^2Ef^-\-Ep—Ey ^Ej^t -\-Ep) 


k 

2 l,,|2 


ma.x{0,E^i-Ek) J\Ei^,-Ep\ 


rmn 

J\Ey 


d|v| 


[Ep - Ek’Y - |v 


1 2 


X f{Ek:)f{Ep,){l - f{Ek)){l - fiEp)) - fiEk)fiEp)il - fiEk'W - /(^pO) • 

(A.14) 

3. Ea{k)Ea{p) EoL{k')Va{p')'- 


Rn 


Gl 


{‘^^?El Jo 


poo pc 

/ dEk' / 

^0 Jmi 


f*Ek-\-Ep 


dEr, 


X 


k 

2 l „|2 


raa.x(0,Ey—Ei^) JTaax{\Ei;—Ep\,\2Ey—Ep—Et^\) 


(i|q| 


{Ep + Ek) — |q 


1 2 


f{EkJf{Ep,){l - f{Ek)){l - fiEp)) - fiEk)fiEp)il - fiEk'))il - f{Ep,)) . 

(A.15) 


4. Va{k)Ua{p) Va{k')Va{p')-- 


R, 


2Gl 


Q :, S , Q ! 




poo p 

/ dEk’ / 

./ 0 t/ rr 


dEr, 


— E/^z ,Ef^/-\-Ep') 


raax{0,Ef.,-Ek) J\E^,-Ep\ 


d|v| 


{Ep - Ek>)‘^ - 


f{EkJf{Epz){l - f{Ek)){l - f{Ep)) - f{Ek)f{Ep){l - f{EkJ){l - f{Ep,)) 

(A.16) 
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5. Va{k)e (p) i^a{k')e {p'): 


R. 


^2 roo 


/*cxD r 

/ dEk' 

Jo Jn 


dEn 


max(me —-Efc+me) 


X 


2 { 2 -KfEl 

j rf|q|<9(|q| - max(|£^fc - |p||, {E^' - |p'| |))6'(min(£’fc + |p|,^fc' + |p'|) - |q|) 

X { 2 xw ± 1 )^((-E'p + Ekf - |qp - "i-e)^ 

+ J t^lvl^'dvl - ma.x{\Ek - |p'll,-Efc' - |p||))6'(min(£;fc + \p'\,Ek> + |p|) - |v|) 

X AxwiiEp - Ek'f - |vp - mlf 
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(A.17) 


6. Ua{k)e+{p) Ua{k')e+{p'): 
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(A.18) 
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(A.19) 
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(A.20) 


We remind the reader again that |pd is not a free parameter, but is constrained by energy 
conservation. 

The integrals over |q|, |v| and |w| can be evaluated analytically, and it turns out that 
they fall into two different functional forms: 


f , , 2\2 2 2 n X^ 

/ dxia — X ) = a X - ax H-h constant, 

J 3 5 

f , , 2\ 

/ dx[a — X ) = ax —^ + constant. 


The remaining two integrals over Ep and E^' must be performed numerically, although as 
discussed in section 2.3 judicious assumptions about certain distribution functions in the inte¬ 
grand make it possible to pre-evaluate some of the integrals only once, rather than evaluating 
them in real time simultaneously with the numerical solution of the QKEs. 


B Repopulation and damping coefficients in the A/S approximation 

We summarise here the full expressions for the dimensionless repopulation and damping coef¬ 
ficients in the A/S approximation discussed in section 2.2. The quantity A = 27r/ f dTIkkfo(k) 
is a normalisation factor. 

Ca,a = ^ J dUkdIlk'dIlp/dUp dE{kp\k'p')'Y^V‘^[iyaik),Ua{p)\i{k'),i{p')]fo{p)fo{k), 

i 

Ca,s = A j dUkdIlk'dUpfdUp dE{kp\k'p') ^ V‘^[^a{k), j{p)fa{k'), j{p')]fo{p)fo{k), 

Ca,„ = Af dUkdUk'dUp/dUp dE{kp\k'p') ^ V^Waik), j{p)fa{k'), j{p')]fo{p)fo{k), 
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Cafl — ^ j dUkdJlk'dUp'dUp dE[kp\k'p')fa{k) 

x[E K (/c), (p) Ii(A:'),«( pO]/ o(^0/o (pO 

i 

+ X] ''^^[^«(^)>Jb)l^«(^0>j(p0]/o(p)(i -/o(p0) > 

j^l/a,ida 

Ca,i = A j dllkdYlk'dllpidllp 5E{kp\k'p')fQ{k) 

X V^K(A;),Pa(p)N(A;'),*(p')]/o(p)(l - k{p') - k{k')) 

i 

+ '^^Wa{k),j{p)Waik'),j{p')]k{k'){kip')-k{p)) 

j^Vafia 

+ ^ V‘^V'a{k),j{p)ka{k'),j{p')]k{p) , 

j = }AoL,i'a 

Ca,2 = ^ j dUkdUk'dUp/dUp 6E{kp\k'k)kik) 

X V'^Wa{k),j{p)ka{k'),j{p')]ikik')kip) - k{k')k{p) - k{p)k{p'))- 

j = Va,9a 

See equations (2.16) and (2.17) for the implementation of these coefficients in the A/S re¬ 
population and damping terms. 
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